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1.   INTRODUCTION 

Human  beings  are  impelled  to  develop  means  for  automatic  image 
processing  for  at  least  two  major  reasons: 

1.  As  man  expands  into  portions  of  the  universe  formerly 
denied  him,  the  need  for  expanding  his  senses  likewise 
increases. 

2.  Our  need  and  our  ability  to  collect  and  transmit  image  data 
are  increasing  at  a  rapid  ratec 

In  fact,  the  data  collection  rate  is  currently  held  at  a  level  well 
below  capacity  simply  because  of  the  limits  imposed  by  processing.  Thus 
the  extremely  large  volume  of  data  generated-  by  a  multi spectral  source 
will  overburden  the  available  computer  resources.  For  example,  it  has 

been  estimated  that  data  rates  of  the  order  of  10   to  10   bits  per 

14 
year  can  be  expected  in  the  near  future.  A  data  rate  of  10   bits  per 

year  is  equivalent  to  covering  the  surface  of  the  earth  every   18  days  in 
six  spectral  regions  with  a  200  foot  digital  sample  step.  The  digital 
processing  of  this  volume  of  data  will  generate  a  large  computer  require- 
ment. 

The  main  objective  of  this  thesis  is  the  study  of  the  specific 
problems  and  procedures  relating  to  the  LANDSAT  multi spectral  imagery-- 
such  as  geometric  correction  of  the  image,  classification  and  clustering 
procedures. 

1.1  Remote  Sensing  of  Earth's  Resources 

Remote  sensing  of  the  environment  is  rapidly  becoming  a  major 


area  of  research  for  engineers  and  scientists  throughout  the  world. 
Satellites  and  high  altitude  aircraft  provide  ideal  platforms  from 
which  earth  resource  data  may  be  gathered.  The  data  is  gathered  in  the 
form  of  several  spectral  images  of  a  particular  area  of  the  earth  under 
observation.  Each  image  represents  the  spatial  distribution  of  electro- 
magnetic energy  as  seen  through  a  given  spectral  "window."  Information 
about  a  particular  area  is  obtained  through  the  study  of  the  spatial 
and  spectral  variations  of  the  data  for  that  area.  Thus  an  important 
subject  before  the  engineering  and  scientific  community  at  the  present 
time  is  the  processing  of  senses  which  represent  tracts  of  the  earth's 
surface  as  viewed  from  above. 

A  typical  scene  consists  primarily  of  regular  and/or  irregular 
regions  arranged  in  a  patchwork  manner  and  each  containing  one  class 
of  surface  cover  type.  These  homogeneous  regions  are  the  "objects" 
in  the  scene.  A  basic  processing  goal  is  to  locate  and  classify  the 
objects  and  produce  a  description  of  the  scene  in  terms  of  tabulated 
results  and/or  a  "type  map." 

As  in  the  image  processing  applications,  the  locations  and 
spatial  features  (e.g.,  size,  shape,  orientation)  of  objects  are  re- 
vealed by  changes  in  average  spectral  properties  that  occur  at  boundaries 
But  unlike  most  other  applications,  the  spatial  features  of  an  object 
often  have  only  a  weak  relationship  to  its  class.  Research  has  shown, 
however,  that  many  classes  can  be  distinguished  reasonably  well  on  the 
basis  of  their  spectral  features,  using  statistical  pattern  classifica- 
tion techniques. 


Remote  sensing  combined  with  a  variety  of  surface,  airborne, 
and  space  platforms  can  now  give  us  an  integrated  overview  previously 
difficult,  if  not  impossible,  to  obtain.  In  turn  we  are  learning  to 
extract  from  these  data  a  new  class  of  information  about  the  earth-- 
information  necessary  for  making  decision  in  what  we  might  term  "planetary 
engineering";  we  foresee  the  development  of  interrelated  physical  methods 
of  the  earth  and  its  environment  that  will  predict  both  the  natural 
course  of  resources  and  environmental  conditions  as  well  as  the  effects 
of  human  actions. 


2.  THE  LANDSAT  PROJECT 

The  primary  objective  of  the  project  [1]  is  to  conduct  experi- 
ments in  the  interpretation  and  utilization  of  earth-resources  survey 
data  from  space.  Doing  this  means  collecting  substantial  quantities  of 
data  about  the  earth's  surface  from  orbital  altitudes,  and  then  trans- 
formitting  it  in  suitable  format  to  the  scientists  who  will  analyze  it 
and  determine  if  it  can  yield  useful  information.  Thus  the  knowledge 
gained  from  the  application  of  data  acquired  by  the  LANDSAT  will  point 
the  way  toward  development  of  fully  operational  and  more  effective  sys- 
tems for  earth  resources  management. 

In  late  July  1972,  the  National  Aeronautics  and  Space  Admin- 
istration (NASA)  launched  a  one-ton  satellite  into  a  sun-synchronous 
orbit  some  560  miles  up  with  a  period  of  103  minutes.  Every  period  the 
satellite  crosses  the  equator  at  about  9:30  a.m.  local  time  on  the  day- 
light side.  The  spin  of  the  earth  during  each  satellite  period  causes 
the  ground  under  the  orbit  to  slip  westward  about  1,550  nautical  miles 
each  period.  After  a  full  day  the  orbit  precision  combined  with  the 
fact  that  the  orbital  period  is  not  an  exact  submultiple  of  the  day, 
insures  that  the  satellite  traces  out  a  track  on  the  ground  in  which  each 
orbital  sweep  lies  about  85  nautical  miles  west  of  the  corresponding 
sweep  of  the  day  before.  Thus  every  18  days,  the  satellite  will  revisit 
nearly  the  same  ground  points,  always  by  morning.  The  meteorological 
cloud-cover  photographs  made  by  the  satellite  are  higher-grade  versions 
of  the  ones  that  are  now  familiar  on  the  daily  television  weather  pro- 
grams; they  have  a  resolution  of  about  one  kilometer  on  the  surface  by 


day  and  five  to  ten  kilometers  even  at  night  with  the  de!ep  infrared. 
2.1  LANDSAT  Mission 

To  achieve  its  broad  objectives,  the  mission  of  the  LANDSAT 
provides  for  the  repetitive  acquisition  of  high  resolution  data  of  the 
earth's  surface  on  a  global  basis.  LANDSAT  carries  three  5+inch  lenses 
that  produce  television  images  of  the  same  scene  in  three  distinct  color 
bands.  This  is  exceptionally  good  television,  with  4,000  lines,  about 
10  times  the  detail  of  home  or  weather-satellite  television.  It  is 
the  RCA  developed  system,  the  Return  Beam  Vidicon  (RBV)  system.  The 
satellite  also  has  a  second  image  system--!' t  is  the  Hughes  Aircraft  de- 
velopment called  Multispectral  Scanner  (MSS).  Here  the  field  is  scanned 
not  electronically  but  mechanically  by  an  oscillating  mirror.  The 
scheme  allows  the  use  of  many  detectors,  one  for  each  of  four  wave  bands, 
without  any  problem  of  registration  since  eyery   detector  samples  one  and 
the  same  scanning  spot.  In  addition,  the  LANDSAT  observatory  will  be 
utilized  as  a  relay  system  to  gather  data  from  remote,  widely  distri- 
buted, earth-based  sensor  platforms  equipped  by  individual  investigators. 
The  data  acquired  by  the  total  LANDSAT  systems  will  thus  permit  quanti- 
tative measurements  to  be  made  of  earth-surface  characteristics  on  a 
spectral,  spatial,  and  temporal  basis. 

The  overall  system  is  illustrated  in  Figure  1.  The  Operations 
Control  Center  (OCC)  is  the  focal  point  of  the  mission  orbital  operations, 
There  the  overall  system  is  scheduled,  spacecraft  commands  are  originated 
and  orbital  operations  are  monitored  and  evaluated.  Data  Collection 
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Figure  1   Overall  LANDSAT  System 


System  (DCS),  telemetry  and  command  data  transfer  between  the  OCC  and 
remote  sites  is  accomplished  by  NASA  communications.  The  NASA  Data  Pro- 
cessing Facility  (NDPF)  accepts  payload  video  data  in  the  form  of  mag- 
netic tapes  received  in  real  time  at  the  NTTF  station  via  the  OCC  or 
by  mail  from  Alaska  and  Goldstone.  The  NDPF  then  performs  the  video-to- 
film  conversion  and  correction,  producing  black  and  white  images  from 
individual  spectral  bands  and  color  composites  from  several  spectral 
bands. 

2.1.1  Observatory  System 

The  elements  of  the  observatory  system  include  the  payload 
systems  and  the  various  support  subsystems  comprising  the  spacecraft 
vehicle.  Control  of  observatory  attitude  to  the  local,  vertical  and 
orbit  velocity  vectors  within  0.7  degree  of  each  axis  is  achieved  by 
a  three-axis  active  attitude  control  subsystem.  It  uses  horizon  scan- 
ners for  pitch  [2]  (pitch  variations  cause  a  change  in  the  vertical 
scale  by  changing  the  vertical  size  of  the  frame.  The  magnitude  is 
Ay  =  he  while  9  is  the  pitch  variation  from  the  top  to  bottom  of  the 
frame)  and  rol 1  (roll  causes  a  skew  in  the  horizontal  direction  of  magni- 
tude Ax  =  he  where  h  is  the  nominal  altitude,  6  is  the  roll  variation 
r  r 

over  the  frame,  and  Ax  is  the  amount  of  horizontal  skew)  control,  and 
a  gyro-compassing  for  yaw  (yaw  variations  are  due  to  errors  in  the  alti- 
tude control  system)  orientation.  An  independent  passive  Attitude  Mea- 
suring Subsystem  (AMS),  operating  over  a  narrow  range  of  about  2  degrees 
provides  pitch  and  roll  attitude  data  accurate  to  within  0.07  degree  to 


8 


aid  in  image  location.     Orbit-adjustment  capability  is  furnished  by  a 
monopropellant  hydrazine  subsystem  employing  one  pound  force  thrusters. 
This  system  is  used  to  remove  launch  vehicle  injection  errors,  and  to 
provide  periodic  trim  to  maintain  a  precise  orbit. 

Payload  video  data  are  transmitted  to  ground  stations  over 
two  wideband  S-band  data  links. 

2.1.2     Payloads 

2 . 1 . 2 o 1  Return  Beam  Vidicon  Camera 

The  RBV  camera  system  operates  by  shuttering  three  independent 
cameras  simultaneously,  each  sensing  a  different  spectral  band  in  the 
range  of  0.48  to  0.83  micrometers.  Since  these  are  visible  wavelengths, 
the  RBV  is  operated  only  in  daylight.  The  viewed  ground  scene,  100  by 
100  nautical  miles  in  area,  is  stored  on  the  photosensitive  surface  of 
the  camera  tube  and,  after  shuttering,  the  image  is  scanned  by  an  elec- 
tron beam  to  produce  a  video  signal  output.  Each  camera  is  read  out 
sequentially,  requiring  about  3.5  seconds  for  each  of  the  three  spectral 
images.  To  produce  overlapping  images  along  the  direction  of  spacecraft 
motion,  the  cameras  are  reshuttered  eyery   25  seconds.  The  video  band- 
width during  readout  is  3.5  MHz.  Orientation  of  the  three  camera  heads 
is  shown  in  Figure  2. 

2.1.2.2  Multi spectral  Scanner 

The  Multispectral  Scanner  (MSS)  is  a  line  scanning  device  which 
uses  an  oscillating  mirror  to  continuously  scan  perpendicular  to  the 
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Camera  Mounted  in 
Spacecraft 


185  km  x  185  km 
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Direction  of 
Flight 


Figure  2  RBV  Camera  Head  Orientation 
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spacecraft  velocity  as  shown  in  Figure  3.  Six  lines,  with  the  same 
bandpass,  are  scanned  simultaneously  in  each  of  the  four  spectral  bands 
for  each  mirror  sweep.  Spacecraft  motion  provides  the  along-track  pro- 
gression of  the  six  scanning  lines.  Optical  energy  is  sensed  simultan- 
eously by  an  array  of  detectors  in  four  visible  spectral  bands  from  0.5 
to  1.1  micrometer  for  daylight  operations  of  LANDSAT.  A  fifth  band  in 
the  near  (thermal)  infrared  from  10.4  to  12.6  micrometers  is  installed. 
The  detector  outputs  are  sampled,  encoded  to  six  bits  and  formatted  into 
a  continuous  data  stream  of  15  megabits  per  second.  During  image  data 
processing  in  the  ground  data  handling  system  facility,  the  continuous 
strip  imagery  is  transformed  to  framed  images  with  a  10  percent  overlap 
of  consecutive  frames  and  an  area  coverage  approximately  equal  to  that 
of  the  RBV  images. 

2.1.3  Wideband  Video  Tape  Recorders 

The  use  of  data  from  the  RBV  and  MSS  sensors  are  complementary 
in  several  respects,  and  both  sensors  are  generally  operated  simultan- 
eously over  the  same  terrain  during  daylight  hours.  When  operated  over 
a  ground  receiving  station,  their  data  are  transmitted  in  real  time  to 
the  ground  receiving  site  and  recorded  there  on  magnetic  tape. 

2.1.4  Data  Collection  System 

The  Data  Collection  System  (DCS)  obtains  data  from  remote, 
automatic  data  collection  platforms,  which  are  equipped  by  specific  in- 
vestigators, and  relays  the  data  to  ground  stations  whenever  the  LANDSAT 
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spacecraft  can  mutually  view  any  platform  and  any  one  of  the  ground  sta- 
tions. Each  DCS  platform  collects  data  from  as  many  as  eight  sensors, 
supplied  by  the  cognizant  investigator,  sampling  such  local  environ- 
mental conditions  as  temperature,  stream  flow,  snow  depth,  or  soil  mois- 
ture. Data  from  any  platform  is  available  to  investigators  within  24 
hours  from  the  time  the  sensor  measurements  are  relayed  by  the  spacecraft, 

2.1.5  LANDSAT  Multi spectral  Scanner  Data 

LANDSAT  multi spectral  scanner  data  is  received  from  the  satel- 
lite by  NASA,  processed,  and  delivered  to  users  recorded  on  computer 
compatible  tape  and  in  photographic  form.  The  computer  tape  form  of 
the  data  is  calibrated  and  line  length  adjusted  by  NASA  but  no  geometric 
corrections  are  applied.  The  system-corrected  photographic  products 
are  corrected  for  many  geometric  distortions  including  earth  rotation 
effects  in  addition  to  the  above  two  corrections.  Also,  these  images 
are  rescaled  so  that  the  horizontal  and  vertical  scales  are  the  same. 
Thus,  the  digital  form  of  the  MSS  data  contains  many  geometric  distor- 
tions and  users  of  this  data  are   faced  with  the  problem  of  compensating 
for  these  errors. 

When  digital  MSS  data  is  reproduced  in  image  form  on  a  stan- 
dard IBM  computer  line  printer  the  resulting  scale  factor  is  approxi- 
mately 1  in.  =  22,400  in.  in  the  horizontal  direction  and  1  in.  = 
24,200  in.  in  the  vertical  direction.  This  scale  differential  exists 
in  addition  to  the  skew  due  to  earth's  rotation  and  all  other  geometric 
errors.  Similarly  when  this  data  is  reproduced  on  a  video  display 
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device  the  horizontal  scale  is  56  meters  per  point  and  the  vertical 
scale  is  79  meters  per  point. 

2.1.6  MSS  Digital  Data  Geometric  Characteristics 

The  LANDSAT  MSS  System  produces  four  spectral  band  digitized 
imagery  of  approximately  100  nautical  mile  (185  km)  wide  strips  beneath 
the  satellite  path.  The  scanner  has  an  instantaneous  resolution  of  79 
meters  and  scan  lines  are  sampled  at  a  rate  such  that  samples  are  spaced 
approximately  56  meters  apart  and  successive  scan  lines  are  spaced  ap- 
proximately 79  meters  apart  as  determined  by  the  forward  motion  of  the 
satellite.  The  image  data  are  edited  so  that  the  along  size  is  approxi- 
mately 96.3  nautical  miles  (155  km).  The  resulting  data  set  consists 
of  nominally  3240  samples  horizontally  (E-W)  and  2340  samples  vertically 
(N-S).  The  geometric  distortions  are  due  to  sensor,  satellite,  and 
earth  effects. 
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3.  COMPILATION  AND  GEOMETRIC  REGISTRATION 
OF  LANDSAT  MULTI TEMPORAL  IMAGERY 


Within  the  context  of  computer-assisted  earth-resources  data 
analysis,  image  registration  objectives  fall  generally  into  two  separate 
(but  procedurally  related)  categories: 

1.  A  first  set  of  image  registration  objectives  has  to  do 
with  the  need  for  accurate  spatial  registration  of  imagery  with  respect 
to  conventional  geographic  coordinate  systems.  Generally,  the  goal  of 
procedures  here  is  to  relate  earth-resources  imagery  of  interest  to  map 
bases  in  such  a  manner  that  they  can  be  associated  with  each  data  ele- 
ment of  a  multispectral  scanner  (MSS)  image  (multivariate  observation, 
resolution  element,  pixel)  an  accurate  ground  position  in  terms  of  a 
conventional  map  coordinate  system,  typically  latitude-longitude.  Such 
image-to-map  registration  is  necessary,  for  example,  to  enable  conven- 
ient cross  referencing  of  image  data  with  other  geographically-referenced 
ground-collected  information  within  automated-classifier  training  opera- 
tions and  other  image  analysis  procedures.  Image-to-map  registration  is 
also  required  where  image  data  are  to  be  output  in  graphical  format 
suitable  for  photo-reproduction  as  topo  map  overlays.  In  this  instance, 
not  only  must  the  mathematical  function  relating  image  coordinates  to 
map  coordinates  be  determined,  but  also  the  data  of  the  image  must  be 
physically  reformatted  such  that  a  graphical  display  of  the  output  data 
via  some  photo-recording  device  conforms  in  area,  orientation,  and  pro- 
jection to  a  specific  map. 

2.  A  second  category  of  objectives  within  image  registration 
operations  relates  to  the  automatic  overlay  of  multiple  images  (usually 
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two)  taken  at  different  times  over  the  same  ground  area  in  creating 
multitemporal  multispectral  image  data  sets.  In  the  case  of  imagery 
collected  by  LANDSAT,  conventional  mathematical  correlation  techniques 
may  be  employed  directly  to  determine  points  of  optimal  match  between 
temporally  separate  images.  The  goal  of  digital  processing  here  of 
course  is  to  determine  automatically  an  element-by-element  matching 
of  the  pixels  of  one  image  with  those  of  the  other  image(s)  such  that 
a  composite  multitemporal  data  set  might  be  produced  for  a  specified 
ground  area  common  to  both  (all)  images. 

Usually,  it  is  also  desirable  that  the  resulting  mutli temporal 
imagery  be  accurately  registered  to  some  specified  geographic  coordinate 
system.  Where  one  map-registered  image  already  exists  for  a  specific 
ground  area,  subsequent  images  can  be  registered  automatically  to  the 
same  map  base  by  simply  correlating  later  images  with  the  initial  map 
registered  image  and  then  applying  identical  image-to-map  spatial  trans- 
formation. Such  combined  uses  of  image- to- image  and  image-to-map  regis- 
tration procedures  suggest  the  practicality  of  automated  systems  for 
monitoring  natural  resource  and  land  use  boundary  changes  over  time. 

The  remainder  of  this  chapter  describes  the  approach  taken  at 
the  Center  of  Advanced  Computation,  University  of  Illinois  [33  in 
adapting  image  registration  procedures  previously  researched  at  LARS 
[2,4,5]  Purdue  University,  for  more  economical  processing  using  the 
ILLIAC  IV  and  other  computers  of  the  ARPA  network. 

3.1  A  Basic  LANDSAT  Data  Preprocessing  Problem:  Image  Rectification 

For  a  number  of  reasons,  it  is  generally  most  convenient  to 
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work  with  the  digital  data  of  LANDSAT  MSS  imagery  in  units  of  map- 
oriented  rectangular  image  areas.  For  example,  the  image  analysis  area 
of  interest  is  often  the  geographic  region  of  a  particular  map  bounded 
by  lines  of  latitude  and  longitude.  However,  the  characteristics  of  the 
near-polar  orbit  of  LANDSAT  are  such  that  the  grid  of  data  elements 
scanned  for  each  MSS  image  is  always  rotational ly  displaced  to  a  consid- 
erable degree  from  alignment  with  geodetic  meridians  and  parallels.  Also, 
the  rotation  of  the  earth  under  LANDSAT  during  the  scanning  of  each  MSS 
image  introduces  an  oblique  skew  to  the  grid  geometry  of  each  image 
digitalized.  Hence,  to  enable  convenient  retrieval,  analysis,  and  dis- 
play of  selected  MSS  imagery  in  terms  of  map-oriented  rectangular  areas, 
cumbersome  digital  preprocessing  is  usually  required  to  transform  geo- 
metrically and  reformat  the  computer  tape  data  of  each  LANDSAT  image  to 
north-vertical  image  grid  orientation. 

This  computer  preprocessing  of  LANDSAT  MSS  data  from  initial 
scanner  format  to  map-oriented  image  analysis  format  is  known  generally 
as  image  rectification  or  geometric  correction.  Again,  the  two  principle 
objectives  of  this  data  preprocessing  step  are: 

1.  To  remove  the  oblique  skew  within  the  MSS  image  grid 
caused  by  the  earth's  rotation  during  image  scanning, 

2.  To  transform  the  MSS  image  grid  to  a  vertical  north,  map 
orientation. 

A  third  objective  to  be  served  by  geometric  correction  is  sometimes 
the  differential  rescaling  of  the  MSS  image  grid  along  its  transformed 
axes  to  achieve  equality  of  horizontal  and  vertical  image  scales. 
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In  any  case  via  image-rectification  data  preprocessing,  a  new 
image  data  tape  is  created  from  which  map-oriented  rectangular  areas  of 
data  can  be  retrieved  and  displayed  directly  for  all  subsequent  image 
interpretation  purposes. 

To  understand  the  nature  of  the  data  processing  problem  as- 
sociated with  geometric  rectification  of  LANDSAT  MSS  imagery,  consider 
the  following  typical  image  rectification  task.  We  wish  to  extract  from 
a  single  LANDSAT  MSS  image  those  data  elements  corresponding  to  ground 
points  lying  within  a  specific  USGS  1°  x  2°  map  quadrangle,  and  to  trans- 
form that  image  subarea  such  that  ordered  rectangular  blocks  of  data  of 
the  rectified  image  tape  conform  closely  in  geometric  area  and  orienta- 
tion for  the  USGS  7-1/2  minutes,  15  minutes,  and  1°  x  2°  topo  maps  exist- 
ing for  that  quadrangle. 

Assume  that  a  mathematical  function  is  given  that,  for  each 
element  of  the  rectified  image  designated  by  output  image  row  and  column 
coordinates  (r  ,  c  ),  determines  that  single  element  of  the  original 
scanner  image  (r,  c)  that  is  most  nearly  geographically  equivalent.  We 
assume  here  an  integer  function,  say  (r,c)  =  f(r  ,c  ),  that  establishes 
an  integer-pair  mapping  from  rectified  pixel  indices  back  to  "best- 
corresponding"  scanner  image  indices.  (Note  that  such  a  mapping  between 
rectified  and  original  image  pixels  will  not  in  general  be  one-to-one 
[next  paragraph  describes  an  oblique  transformation  method  to  achieve 
a  one-to-one  mapping  of  LANDSAT  input  picture  grid  onto  a  LANDSAT  output 
picture  grid]  and  thus  some  data  elements  are  usually  lost  or  duplicated 
in  tran format ion.) 


18 


Having  a  particular  rule  for  locating  pixels  of  the  original 
image  to  be  copied  into  specific  pixel  positions  within  the  rectified 
image  format,  we  then  wish  to  apply  this  reformatting  rule  for  all  pixels 
lying  within  the  given  1°  x  2°  map  quadrangle.  Figure  4  illustrates  the 
impractical ity  of  core-contained  computational  strategies  for  this 
problem. 

We  would  like  to  input  from  magnetic  tape  successive  portions 
of  the  oribinal  scanner  image,  and  following  transformation,  output  se- 
quentially onto  another  tape  individual  complete  lines  of  the  rectified 
image.  The  main  difficulty  arises  with  the  size  of  the  real -core  memory 
required  to  hold  that  portion  of  the  input  image  containing  all  data  ele- 
ments needed  to  compose  the  next  line  of  the  output  image.  Rectifica- 
tion of  the  1°  x  2°  quadrangle  of  LANDSAT  data  as  depicted  in  Figure  4 
would  require  typically  a  real  core  input  data  buffer  on  the  order  of 
7.10  bytes . 

Confronted  with  the  sheer  quantity  of  data  comprising  each 
LANDSAT  MSS  image  and  given  such  computational  problems,  most  LANDSAT 
data  analysts  reexamine  the  necessity  of  preprocessing  LANDSAT  image 
tapes  to  achieve  map-oriented  data  format.  At  this  point,  it  is  usually 
decided  that  the  goals  of  image-to-map  registration  should  be  limited  to 
only  the  determination  of  those  numerical  transformation  necessary  for 
referencing  individual  data  elements  by  map  coordinates  leaving  completely 
undone  the  actual  reformatting  of  LANDSAT  image  data  tapes  to  map  orien- 
tations. In  many  instances  where  geo-referencing  and  graphical  output  is 
only  identical  to  data  analysis  objectives,  such  a  strategy  seems  entirely 
appropriate. 
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Rectified  1°  x  2°  Quad  Image 


Figure  4  Buffer  Requirement  for  Core-Oriented  Rectification  of 
LANDSAT  Images 
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However,  a  fundamental  problem  remains 0  In  order  to  determine 
adequately  (to  one-pixel  average  accuracy)  the  numerical  parameters 
needed  for  accurate  geo-referencing,  some  initial  set  of  control  points 
must  be  identified  simultaneously  in  terms  of  both  LANDSAT  image  pixel 
position  and  map  position.  Such  geographic  control  points  most  conven- 
iently obtained  for  LANDSAT  imagery  over  the  U.S.  by  manually  correlat- 
ing small  area  line-printer  displays  of  rectified  LANDSAT  image  data 
with  USGS  7-1/2  minute  topo  maps  and  digitizing  some  small  set  of  points 
denoting  overlay  positions  of  maximum  geographic  correspondence. 

Here  a  reasonable  approach  might  seem  simply  to  retrieve  from 
tape  and  rectify  only  those  small  areas  of  image  data  needed  for  manual 
correlation  with  quadrangle  maps.  However,  the  rectangular  latitude- 
longitude  grid  of  the  USGS  map  series  is  such  that  a  uniformly  distri- 
buted set  of  image  sub-areas,  selected  to  fall  within  7-1/2  minute  map 
quadrangles,  is  such  more  efficiently  retrieved  from  an  image  data  tape 
of  recti fied-image  format.  Since  also  there  will  be  needed  some  sys- 
tematic grid  of  image  subareas  for  automated  digital  correlation  within 
any  subsequent  mul ti temporal  image  registration  procedures,  it  would 
seem  that  a  number  of  image  registration  objectives  would  be  served  by 
some  efficient  computational  strategy  for  initial  transformation  of  all 
LANDSAT  imagery  to  vertical-north  orientation.  Such  an  efficient  LANDSAT 
data  preprocessing  operation  is  provided  by  the  oblique  (skew)  trans- 
formation. 
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3.2  Oblique  Transformation  of  LANDSAT  Images  to  Approximate  North- 
South  Orientation 


The  oblique  transform  [6]  achieves  a  one-to-one  mapping  of  a 
LANDSAT  input  picture  grid  onto  a  LANDSAT  output  picture  grid  that 
closely  approximates  a  skewed  rotation  of  the  input  picture  with  scale 
factored  out.  A  one-to-one  mapping  is  desirable  since  a  further  cor- 
rection step,  dependent  on  a  particular  map  projection  and  some  set  of 
specific  image-to-map  control  points,  will  be  applied  to  the  obliquely- 
transformed  output  picture.  The  major  purpose  of  the  oblique  transform 
is  to  reduce  the  size  of  the  random-access  memory  needed  for  application 
of  that  final  correction  step.  As  a  bonus,  it  turns  out  that  a  line 
printed  display  (assuming  10  characters  per  inch  horizontally  and  8 
characters  per  inch  vertically)  of  the  obliquely  transformed  LANDSAT 
image  has  a  scale  and  geometry  quite  close  to  that  of  the  USGS  seven- 
and-a-half-minute  quadrange  map  series.  The  oblique  transformation  is 
equivalent  to  a  vertical  skew  of  the  input  picture  followed  by  a  hori- 
zontal skew.  See  Figure  5.  Each  input  column  is  shifted  an  integral 
number  of  pixels,  the  shift  distance  being  the  rounded  off  values  of 
-  3  x  for  that  column.  The  resulting  rows  are  now  shifted  horizontally 
by  the  rounded  off  value  of  ay  for  that  row.  It  is  clear  that  with  this 
process,  no  input  points  are  thrown  away  or  duplicated.  The  question  is 
now  the  relationship  between  the  oblique  transform  and  a  general  linear 
transform. 

The  two  parameters  of  the  oblique  transform  are  the  skew  fac- 
tors a  and  3.  Since  all  input  and  output  pixels  must  fall  within  integral 
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Figure  5  The  Oblique  Transformation  Process 
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grid  location,  the  integral  or  "digitized"  oblique  transform  is: 

x  =  x  -  ay  +  1/2 

'   I   '      I  (1) 

y  =  [Bx  +  1/2J+  y 

with  x  and  y  output  pixel  coordinates,  x  and  y  input  pixel  coordinates, 
and  where  |_aj  denotes  the  greatest  integer  ^  a.  Thus,  [a_  +  1/2J  is  theo- 
retically equivalent  to  rounding  a^  to  the  nearest  integer.  The  trans- 
formation (1)  is  the  inverse  of  the  process  displayed  in  Figure  5  to 
reflect  the  way  picture  transformations  are  normally  obtained:  for  each 
output  pixel,  the  input  pixel  that  transforms  to  it  is  found  and  assessed. 
The  continuous  equivalent  of  (1)  is  the  linear  transform: 

1 
x  =  x  -  ay 

(2) 

y  =  Bx  +  (1  -  a3)  y 

Note  that  rounding  off  x  and  y  values  in  (2)  does  not  necessarily  re- 
sult 1n  an  information-lossless  transform.  For  example,  with  a  =  1/2  and 

1   1 
B  =  -  2,  both  coordinates  (1,1)  and  2,2)  are  transformed  into  (x  ,y  )  = 

(l,o). 

The  general  linear  transform  may  be  factored  into  an  oblique 
transform  (2)  followed  by  a  scale  transformation.  One  factorization  is: 


(3a) 


If  there  is  no  solution  to  (3a),  i.e.,  a  =  0  then  the  alternate  factor- 
ization is: 
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(3b) 


Given  any  non-singular  linear  transform,  i.e.,  ad-bc  f   o,  a,  3>  m  and  n 
ncessarily  exist  and  can  be  easily  determined.  For  more  details  see 
reference  [2], 

3.2.1  Error  Bounds 

The  maximum  errors  incurred  by  using  (1)  instead  of  (2)  are: 

6x  =  ay  -  [ay  +  1/2J  (4a) 

6y   =  3(x  -  ay)  -  [g(x  -  [ay  +  1/2_J  )  +  1/2]  (4b) 

Any  real  number  (u)  can  be  expressed  as  m  +  e,  where  m  is  an  integer, 
and  -1/2  ^  e  <  1/2.  Thus,  we  can  set  ay  =  m  +  e,  and  obtain  from  (4a): 

<5x  =  m+£-[m+£+l/2j   =m+e-m  =  e 
Thus,     J  6x  J  1   1/2.   For  (4b),  we  set  x  -  ay  =  m  +  e-| 
where     -  1/2  <  £]  -  1/2  instead.  Then 

x  -  [ay  +  1/2]   =  x  -  [x-m-E1  +  1/2J  =  m,  and         (5a) 

we  get    6y  =  3m  -  3m  +  1/2   +  3£i  =  e  +  $e, 

1 —     — '       '         '  (5b) 

|5y|  <  1/2  (|3|+D 

From  theoretical  consideration  of  the  LANDSAT  orbit,  the  actual  3  re- 
quired to  bring  LANDSAT  images  over  the  continental  United  States  into 
north-south  alignment  does  not  exceed  .17  in  absolute  value,  giving  a 
maximum  of  error  of  0.59  of  a  pixel's  vertical  dimension. 
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3.2.2  Implementation  Consideration 

Magnetic  tape  used  at  1600  bps  is  currently  the  only  reasonable 
way  to  store  LANDSAT  images.  Small  disk  packs  (IBM  2314  equivalent)  will 
not  hold  a  complete  LANDSAT  picture  and  large  ones  (IBM-  3330  equivalent) 
are  relatively  scarce  and  expensive.  Since  typically  images  are  mailed 
via  tape  anyway,  the  present  implementation  of  the  oblique  transformation 
procedure  assumes  magnetic  tape  storage  of  both  the  output  and  the  input 
images.  The  main  restriction,  therefore,  is  that  all  access  to  the  data 
of  both  input  and  output  pictures  is  necessarily  sequential ,  i.e.,  a 
row  or  line  at  a  time. 

Whenever  the  3  parameter  of  the  oblique  transformation  (1)  is 
non-zero,  the  information  of  more  than  one  input  line  must  be  more  or 
less  randomly  accessible  for  the  creation  of  one  output  line  of  data. 
That  means  input  lines  must  be  buffered.  The  size  of  that  buffer  is 
clearly  the  maximum  number  of  input  lines  needed  for  any  output  line 
multiplied  by  the  memory  needed  per  input  line,  i.e.,  (3.  width)  (width) 
(Kbytes/pixel).  For  a  typical  U.  S.  LANDSAT  frame,  3  =  0.17,  width  = 
3240  pixels  per  row,  and  K  =  4,  8-bit  bytes  per  pixel,  yielding  a 
required  buffer  size  of  slightly  more  than  7.10  bytes.  Note  if  we 
perform  the  transformation  or  quarter-frames  rather  than  the  entire 
image,  all  operations  are  done  four  times,  but  memory  requirements  will 
be  reduced  by  a  factor  of  approximately  16.  Furthermore,  computation 
time  using  such  multiple  passes  will  not  be  significantly  higher  since 
the  same  number  of  inner-loop  operations  are  done  in  either  case. 
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Specific  implementation  of  the  oblique  transformation  pro- 
cedure should  be  designed  with  respect  to  "real  core"  availability,  since 
assessing  the  full  7M  bytes  buffer  as  "virtual  core"  on  most  virtual 
machines  would  create  a  large  number  of  page  transfers  or  faults.  The 
great  increase  in  system  overhead  that  causes  would  reduce  enormously 
the  efficiency  of  the  procedure. 

It  is  interesting  to  compute  the  approximate  number  of  page 
faults  that  would  occur;  each  output  line  causes  at  least  3-  width-w 
faults,  when  w  is  the  working  set  size  or  the  number  of  pages  storable 
in  real  core.  If  we  assume  2340  output  lines,  use  the  same  width  and 
3  as  before,  and  use  w  =  100  (equivalent  to  a  400K  byte  working  set  in 
real  core  on  an  IBM  360/67),  about  1,300,000  faults  would  occur. 

Finally,  even  if  the  input  width  is  cut  down  as  above,  we 
are  still  faced  with  the  condition  that  tape  input  and  output  is  more 
efficient  as  a  real-core  machine,  especially  when  large  tape  input  and 
output  buffers  are  used. 

3.2.3  Details  of  the  Real -Core  Machine  Implementation 

The  process  of  oblique  transformation  is  broken  up  into  several 
processing  steps,  each  step  transforming  a  successive  part  of  a  LANDSAT 
input  picture.  Thus,  each  pass  has  one  output  (the  freshly-updated  out- 
put picture)  and  two  inputs  (the  appropriate  LANDSAT  MSS  quarter-frame 
tape  as  received  from  Goddard  and  the  output  picture  from  the  immediately 
preceding  pass). 
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The  part  of  a  LANDSAT  tape  that  is  transformed  in  each  pass  is 
a  vertical  strip  of  it,  e.g.,  the  entire  Goddard  tape  if  sufficient  real 
core  is  available.  Within  each  step  then,  for  each  output  line,  the 
following  things  are  done.  First,  the  output  line  is  "clipped"  into  the 
largest  segment  that,  transformed,  lies  wholly  within  the  current  input 
strip. 

Second,  if  that  segment  is  non-empty,  then  the  required  lines 
of  the  input  picture  are  read  in  and  queued  into  the  "back"  part  of  the  in- 
put buffer.  Note  that  due  to  the  linear  nature  of  the  oblique  transfor- 
mation, lines  previously  queued  which  are  not  needed  for  the  current  output 
line  will  necessarily  be  at  the  "front"  part  of  the  input  buffer  and  can 
be  discarded  systematically  with  the  knowledge  that  their  information 
will  not  be  needed  again  in  the  current  pass.  That  fact  allows  us  to 
efficiently  keep  the  size  of  the  queue  (the  input  lines  buffer  as  dis- 
cussed above)  just  that  needed  to  hold  all  input  line  data  for  a  single 
output  segment. 

Finally,  the  accumulated  data  for  the  current  output  line  is 
read  from  the  tape  output  from  the  last  pass,  and  for  each  point  in  the 
output  segment,  the  inner  loop  of  the  oblique  transformation  updates 
that  part  of  the  output  line  by  moving  the  corresponding  image  data  from 
the  input  line.  The  new  output  line  is  then  written  out  onto  another 
tape. 

3.2.4  Precision  Geographic  Registration  of  LANDSAT  Imagery 

Initial  rectification  of  all  LANDSAT  imagery  to  map-oriented 
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data  formats  greatly  facilitates  all  subsequent  precision  image-to-map 
registration  efforts.  Once  approximate  registration  has  been  determined 
for  an  image  and  data  tapes  re-formulated  via  skew  transformation,  pre- 
cision image-to-map  registration  may  be  conveniently  accomplished  inthe 
following  manner: 

The  central  problem  confronting  all  precision  image-to-map 
registration  efforts  is  to  identify  accurately,  in  terms  of  both  digital 
image  pixel  coordinates  and  map  coordinates,  a  geographic  network  of 
contral  points  visible  both  on  the  image  and  on  maps.  For  LANDSAT  image- 
ry over  the  U.S.,  experience  suggests  that  such  a  system  of  control 
points  can  be  established  by  overlaying  on  a  light  table  small  line- 
printer  displays  of  rectified  imagery  with  USGS  7  1/2  quadrangle  maps 
at  1:24,000  to  determine  manually  image-to-map  overlay  positions  of 
maximum  geographic  correspondence. 

The  fact  that,  the  horizontal  and  vertical  scales  of  standard 
line-printer  displays  of  rectified  LANDSAT  imagery  (8  pixels  to  the 
inch  vertically  and  10  pixels  per  inch  horizontally)  overlay  USGS  7  1/2 
minute  quadrangles  to  approximately  one-pixel  accuracy  contributes  much 
toward  the  convenience  of  such  a  strategy.  Following  quick  visual  cor- 
relation of  line  printer  display  and  map  quadrangle,  a  single  point  of 
geographic  correspondence  between  image  and  map  can  be  established  simply 
by  recording  the  row  and  column  coordinates  of  a  central  pixel  of  the 
display  and  measuring  the  position  of  the  center  of  this  pixel  in  terms 
of  map  coordinates.  The  availability  of  a  data-tablet  digitizer  for 
such  a  procedure  greatly  facilitates  the  establishment  of  a  network  of 
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control  points  by  this  method.  In  selecting  image  areas  for  determina- 
tion of  control  points,  essentially  two  strategies  are  available. 

We  have  noted  that  the  skew-transformed  image  tape  now  stores 
approximate  map  quadrangles  of  imagery  as  ordered  rectangular  blocks  of 
data  since  we  will  need  some  grid  of  image  blocks  (64  lines  by  64  col- 
umns) for  any  subsequent  image-to-image  registration  procedures  that  may 
be  required.  It  is  convenient  in  several  respects  to  establish  for  each 
rectified  image  a  grid  of  image  correlation  blocks  centered  approximately 
on  USGS  7  1/2  minute  quadrangles  to  be  used  throughout  all  image-to-map 
and  image-to-image  registration  procedures.  Once  this  quad-centered 
grid  of  image  blocks  has  been  determined  and  retrieved  from  tape,  line- 
printer  displays  (MSS  channel  two)  can  be  produced  for  only  a  small  addi- 
tional cost.  Quick  visual  inspection  of  these  displays  with  reference 
to  available  quadrangle  maps,  should  suggest  some  subset  of  image  blocks 
to  be  used  in  establishing  precision  geographic  control.  However,  such 
a  strategy  is  appropriate  only  for  those  geographic  regions  for  which 
USCS  7  1/2  mapping  is  relatively  complete. 

In  those  geographic  regions  for  which  7  1/2  quadrangle  maps 
are  unavailable,  image-to-map  registration  must  proceed  in  more  conven- 
tional fashion.  Some  sets  of  candidate  control  points  must  be  selected 
manually  by  the  data  analyst  with  reference  to  suitable  terrain  features 
identifiable  on  both  film  image  and  available  maps.  A  set  of  small 
image  blocks  containing  these  candidate  control  points  is  retrieved  from 
tape  and  displayed  in  line-printer  image  format.  Then,  whenever  it  is 
meaningful  to  do  so  given  available  materials,  individual  image  pixels 
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are  recorded  in  terms  of  both  yow  and  column  positions  and  geographic 
coordinates. 

Our  initial  experience  suggests  that  to  achieve  one-pixel 
average  registration  accuracy  over  a  full  LANDSAT  image,  approximately 
two  dozen  widely  distributed  control  points  should  be  established  by 
either  of  those  methods. 

3.2.5  Compilation  of  LANDSAT  Multitemporal  Imagery 

In  overlaying  multiple  LANDSAT  images  to  achieve  multitemporal 
data  sets,  points  of  best  match  between  the  two  images  to  be  overlaid 
must  be  determined.  These  match  points  may  be  determined  digitally  by 
computing  points  of  maximum  image  correlation  between  small  blocks  of 
data  paired  between  the  two  images.  Following  the  LARS  image- to- image 
registration  methodology  [5]  such  a  procedure  may  be  implemented  as  fol- 
lows. 

Given  two  temporally  separate  LANDSAT  images  over  the  same 
ground  area,  initial  geographic  control  allows  pairs  of  image  correla- 
tion blocks  to  be  selected  from  the  two  images  such  that  spatial  dis- 
placement between  paired  blocks  are  only  on  the  order  of  10  pixels. 
One  "floater"  block  of  data  taken  from  one  image  is  translated  hori- 
zontally and  vertically  through  all  possible  positions  of  pixel-to- 
pixel  overlay  with  the  corresponding  (and  slightly  larger)  "stationary" 
block  taken  from  the  other  image.  Numerical  values  of  image  correla- 
tion are   computed  for  all  overlay  positions,  and  thus  a  two-dimensional 
correlation  surface  is  determined.  A  single  point  of  spatial  correspond- 
ence between  the  two  image  blocks  is  then  established  automatically  at 
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the  peak  of  this  surface.  Repeated  for  a  complete  grid  of  blocks  paired 
between  two  images  such  a  procedure  outputs  a  computer-determined  system- 
atic set  of  control  points  which  can  then  be  employed  directly  in  cali- 
brating the  spatial  transformation  required  for  overlaying  the  two  images 
as  a  single  multi temporal  image. 

For  typical  applications,  it  is  sufficient  to  choose  floater 
blocks  from  one  image  as  32-line  by  32-column  data  rectangles  and 
stationary  blocks  from  the  other  images  as  64  by  64  rectangles.  Allow- 
ing all  possible  translations  of  floaters  with  respect  to  stationary 
blocks,  1089  (33  x  33)  individual  image  correlations  are  computed  for 
each  block  pair. 

Also,  it  is  sufficient  to  space  correlation  blocks  approximate- 
ly 180  pixel  positions  apart,  both  horizontally  and  vertically,  i.e., 
the  approximate  dimensions  of  a  7  1/2  minute  map  quadrangle.  Therefore, 
for  two  LANDSAT  images  perfectly  coincident  over  the  same  geographic 
area,  correlation  would  be  computed  for  approximately  230  block  pairs. 
This  corresponds  to  more  than  250,000  individual  image  correlation  cal- 
culations. Since  there  image  correlation  calculations  represent  the 
majority  of  all  computations  needed  for  registration  of  two  images,  and 
since  these  calculations  can  be  performed  in  parallel  computation  for- 
mat, it  is  convenient  to  isolate  these  correlation  calculations  for 
execution  on  the  ILL (AC  IV.   This  can  be  done  by  retrieving  from  both 
image  tapes  of  UCLA  two  quad-centered  grids  of  image  correlation  blocks 
paired  between  the  two  images  and  transmitting  these  two  files  via  the 
ARPA  network  to  ILLIAC  IV  at  NASA/ Ames.  Following  the  ILLIAC  IV  calcu- 
lations of  all  block  correlation  peaks,  only  minor  TENEX  processing  is 
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required  to  calculate  the  numerical  parameters  of  the  two-dimensional 
least  squares  quadratic  polynomial  needed  for  pixel-to-pixel  overlay  of 
the  images.  Given  the  parameters  of  this  polynomial,  the  UCLA  IBM  360/91 
may  be  again  employed  efficiently  for  re-formatting  the  two  separate 
images  as  a  single  8-channel  multi temporal  image. 

3.2.6  The  EDITOR.  .  .ERTS  Data  interpreter  and  TENEX  Operation  Recorder 

The  EDITOR  system  represents  a  component  within  a  more  exten- 
sive computer  software  development  effort  undertaken  at  CAC  to  facili- 
tate large-scale  ILLIAC  IV  -  ARPA  Network  analysis  of  multispectral 
reconnaissance  imagery  such  as  that  collected  by  the  Earth  Resources 
Technology  Satellite--ERTS  (currently  called  LANDSAT)  of  NASA  and  USGS/ 
DI.  [9,10]  This  software  effort  has  been  undertaken  by  CAC,  in  collab- 
oration with  the  Laboratory  for  Applications  of  Remote  Sending  (LARS) 
of  Purdue  University,  in  support  of  the  LANDSAT/EROS  programs  of  NASA 
and  USGS/DI,  in  support  of  the  ILLIAC  IV  and  ARPA  Network  applications 
programs  of  NASA  and  ARPA/DOD,  and  in  support  of  U.  S.  agricultural 
crop-acreage  monitoring  objectives  of  SRS/USDA.  Image  analysis  algorithms 
initially  implemented  within  these  efforts  follow  closely  software  pre- 
viously researched  at  LARS  as  multispectral  image  interpretation  tech- 
niques. [11,12] 

EDITOR  has  been  conceived  as  an  interactive  job-editing  and 
file-management  interface  to  ERTS  data  analysis  systems  being  implemented 
for  the  ILLIAC  IV  hardware  complex  at  NASA's  Ames  Research  Center  (Ames) 
at  Moffett  Field,  California.  Nationally  accessible  via  the  ARPA  Network, 
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the  EDITOR  system  allows  portable,  dial-up  terminal  command  of  ILLIAC 
IV  image  processing  facilities  at  Ames.  Multispectral  image  data  manage- 
ment services  within  EDITOR  have  been  designed  to  optimize  the  efficiency 
of  ERTS  data  file  transfers  between  the  computers  of  the  ILLIAC  IV  com- 
plex at  Ames  (the  ILLIAC  IV  parallel  processor,  its  peripheral  TENEX 
PROCESSORS,  AND  THE  UNICON  Data  Computer)  and  other  computers  on  the  ARPA 
Network. 

Since  the  present  EDITOR  system,  using  tape  inputs,  has  been 

developed  for  stimulation  of  data  management  systems  to  be  implemented 

12 
using  the  10  -bit  laser  memory  of  the  UNICON  Data  Computer  at  Ames, 

specific  data  tape  formatting  is  required  for  all  multispectral  images 

to  be  processed  by  EDITOR.  Required  tape  formats  are   described  in  CAC 

Technical  Memorandum  No.  19.  [13] 

While  conceived  initially  as  a  peripheral  TENEX  file-management 
and  job-editing  interface  to  ILLIAC  IV  batch  image  interpretations,  the 
scope  of  EDITOR  functions  has  been  widened  subsequently  to  include  pro- 
cedures for  direct  TENEX  interpretation  of  small-scale  ERTS  data  samples 
in  preparation  for  ILLIAC  IV  batch  analyses.  Thus,  where  only  small- 
scale  image  interpretations  are  needed  and  where  portable  terminal  output 
is  sufficient,  EDITOR  may  be  used  alone  as  a  self-sufficient  tape- 
operating  multispectral  image  analysis  system. 

Since  EDITOR  has  been  developed  as  an  interactive  interface  to 
ILLIAC  IV  image  processing  facilities,  the  system  should  also  prove 
advantageous  to  a  wider  community  of  ARPA  Network  picture  processing 
researchers.  Where  image  disk  files  conform  to  EDITOR  data  file  formats, 
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the  system  may  be  used  for  image  file  management,  editing,  and  display 
in  conjunction  with  more  general  ILLIAC  IV  image  processing  activities. 
Specific  file  formats  assumed  by  EDITOR  are  described  in  CAC  Technical 
Memorandum  No.  19.  [13] 

The  EDITOR  system  described  here  may  be  executed  on  any  ARPA 
Network  TENEX  system  having  7-track  or  9-track  magnetic  tape  drives. 
Currently,  a  copy  of  the  EDITOR  system  is  maintained  by  CAC  at  Bolt, 
Beranek  and  Newman  (BBN)  in  Cambridge,  Massachusetts,  in  directory 
<RAY>.  The  BBN  EDITOR  can  be  made  available  to  any  user  having  access 
to  the  ARPA  Network.  The  file  <RAY>TAPES. SHARED  is  maintained  there  by 
CAC  describing  all  tape  imagery  available  at  BBN  for  experimental  system 
usage. 

3.2.7  System  Overview 

The  actual  use  of  EDITOR  is  best  described  by  way  of  examples. 
The  following  verbal  description  of  EDITOR  simply  sketches  major  system 
functions.   (More  detailed  descriptions  of  all  presently  available 
system  procedures  are  given  in  Section  3.2.9.)  A  familiarity  with  basic 
multispectral  image  interpretation  procedures  is  assumed. 

EDITOR  is  "implemented  in  modular  format  with  each  system 
module  addressed  by  a  single  command.  For  each  command,  only  that  num- 
ber of  characters  sufficient  to  distinguish  the  command  from  all  others 
need  be  typed  by  the  user. 

The  RETRIEVE  command  is  used  to  create  TENEX  disk  files  con- 
taining multispectral  image  data  from  single  rectangular  areas  of  an 
LANDSAT  image,  or  from  sets  of  rectangular  areas  within  an  image.  Use 
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of  RETRIEVE  assumes  that,  prior  to  entering  EDITOR,  the  user  has  request- 
ed the  TENEX  operator  to  mount  a  specific  image  tape.  Analysis  areas  to 
be  retrieved  from  tape  may  then  be  specified  within  EDITOR  by  the  user 
in  terms  of  image  line  and  column  coordinates  and  line  and  column  sampl- 
ing increments. 

If,  on  the  other  hand,  the  user  has  an  image  file  already 
resident  on  the  TENEX  disk  and  wishes  to  use  only  some  portion  of  it  for 
analysis,  he  may  use  the  CLIP  A  WINDOW  command.  This  command  allows  the 
user  to  construct  new  windows  by  specifying  new  image  line  and  column 
coordinates  for  the  subwindows  desired.  Additional  image  sampling  (by 
additional  line  and  column  sampling  increments)  is  also  possible  within 
CLIP. 

Once  image  analysis  areas  have  been  selected  and  stored  on  the 
TENEX  disk,  the  PRINT  command  within  EDITOR  allows  a  user  to  view  the 
image  disk  files  created  to  verify  that  the  desired  data  windows  have 
indeed  been  retrieved.  Gray-scale  displays  of  any  of  the  individual 
spectral  bands  of  data  within  an  image  file  may  be  achieved  by  overprint- 
ing characters  on  the  user's  terminal.  A  data  histogramming  routine  is 
available  for  displaying  the  proportional  distribution  of  pixel  intensi- 
ties within  any  spectral  band  of  an  image  file.  This  feature  of  the 
system  can  be  used  for  terminal  display  enhancement  purposes.  With 
reference  to  image  data  histograms,  a  user  can  select  manually  a  partic- 
ular range  of  image  point  intensities  to  be  assigned  to  each  gray-scale 
shade  available  within  the  character-overprinting  terminal  display  facil- 
ity. 
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Once  a  user  has  verified  that  the  desired  imagery  has  been 
selected,  several  options  are  available  for  data  analysis. 

The  CLUSTER  command  within  EDITOR  evokes  a  multivariate  cluster 
analysis  of  all  multispectral  image  data  points  within  a  particular  file 
to  determine  a  logical  partitioning  or  grouping  of  all  data  points  into 
a  user-specified  small  number  of  spectral  clusters.  [7,14]  The  cluster 
analysis  procedure  may  be  regarded  as  an  unsupervised  image  interpre- 
tation procedure  categorizing  a  set  of  image  points  into  most-separable 
spectral  classes.  Also,  following  any  cluster  analysis  procedure,  a  set 
of  summary  statistics  characterizing  the  multivariate  distributions  of 
data  within  all  clusters  (i.e.,  mean  vectors  and  variance-covariance 
matrices)  is  computed  and,  at  the  request  of  the  user,  displayed  on  the 
terminal  device. 

The  CLASSIFY  command  within  EDITOR  allows  statistical  classifi- 
cation of  all  data  points  within  a  specified  disk  file  into  a  particular 
set  of  spectral ly-distinct  categories  determined  by  a  previous  cluster 
analysis. 

This  methodology  of  classification  assumes  that  "ground-truth" 
samples  will  always  be  cluster-analyzed  prior  to  class  signature  compu- 
tation in  order  to  verify  the  existence  of  spectral  homogeneity  within 
nominal  terrain  classes  and  spectral  separability  between  classes.  Where 
spectral  properties  of  ground- truth  samples  are  not  homogeneous  within 
nominal  classes  (i.e.,  multivariate  distributions  are  non-spherical),  the 
multiple  spectral  classes  determined  by  cluster  analysis  for  such  terrain 
classes  may  be  used  for  analysis  purposes  as  spectrally-distinct,  but 
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nominally-synonymous,  classes  to  be  re-grouped  following  classification 
into  single  nominal  terrain  classes. 

In  accordance  with  the  LARS  multispectral  image  interpretation 
methodology,  the  statistical  classification  algorithm  evoked  by  the  com- 
mand CLASSIFY  employs  the  Gaussian  maximum-likelihood  decision  rule. 
[11,7,8]  For  each  data  point  of  a  file  of  data  points  to  be  classified, 
discriminant  functions  for  all  classes  are  computed  and  the  point  is 
assigned  to  that  class  for  which  the  discriminant  is  largest.  The  param- 
eters of  the  discriminant  function  for  each  class  are  the  estimates  of 
the  means,  variances  and  covariances  of  the  spectral  properties  of  all 
data  points  belonging  to  that  class.  These  are  the  spectral  class  sta- 
tistics computed  and  saved  by  EDITOR  immediately  following  any  cluster 
analysis  for  the  set  of  spectral  classes  determined  by  that  particular 
CLUSTER  procedure. 

The  ENTER  STATISTICS  FILE  EDITOR  command  provides  a  set  of 
file-management  procedures  convenient  for  interactive  display  and  edit- 
ing of  the  spectral  class  statistics  accumulated  from  previous  cluster 
analyses.  Within  this  mode,  class  signatures  of  different  statistics 
files  may  be  re-combined  in  any  manner  to  produce  new  composite  sets  of 
signatures  to  be  used  for  classification.  Also  within  this  mode,  meas- 
ures of  spectral  separability  or  distinctness  may  be  computed  and  dis- 
played for  all  pairs  of  signatures  within  any  statistics  file.   A  pro- 
cedure is  also  implemented  (again  following  LARS  methods)  for  "pooling" 
the  statistics  of  several  different  signatures  into  a  single  new  compos- 
ite signature. 
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The  EXECUTE  ANALYSIS  ON  ILLIAC  IV  command  instructs  EDITOR  to 
use  the  ILLIAC  IV  at  Ames  for  a  particular  cluster  analysis  or  classifi- 
cation task  specified  by  a  user.  Given  this  command,  EDITOR  immediately 
requires  the  user's  I4-TENEX  user  code  and  password.  If  such  a  code  can 
be  supplied,  EDITOR  will  automatically  set  up  the  task  for  ILLIAC  IV 
processing,  transmit  the  job  via  the  ARPA  Network  to  Ames,  and  enter  the 
job  into  the  ILLIAC  IV  batch  stream.  Procedures  are  also  provided  within 
this  mode  for  inquiring  on  the  status  of  ILLIAC  IV  jobs  and  for  conven- 
ient retrieval  by  EDITOR  of  cluster  analysis  and  classification  output 
files  resulting  from  previous  ILLIAC  IV  executions.  Specific  ILLIAC  IV 
mul tispectral  image  cluster  analysis  and  classification  procedures  acces- 
sible to  EDITOR  users  have  been  documented  elsewhere.  [15,16,17,18] 

In  addition  to  the  procedures  for  displaying  raw  mul tispectral 
imagery,  the  PRINT  command  within  EDITOR  also  provides  procedures  for 
displaying  the  interpreted  imagery  that  results  from  either  cluster 
analysis  or  classification  operations.  Interpreted  image  data  files  may 
be  output  as  terminal  printer  maps  where  the  symbols  1  -  9  and  A  -  Z  are 
used  to  represent  terrain  classes.  Display  of  only  selected  subsets  of 
classes  is  permitted.  Following  again  the  LARS  methodology,  provisions 
are  included  for  reliability  thresholding  of  classification  results  with- 
in PRINT  mode.  Where  a  particular  confidence  level  is  specified  by  the 
user,  EDITOR  displays  only  those  image  data  points  whose  classification 
reliability  meets  the  specified  level.  Also,  procedures  are  included 
within  PRINT  mode  to  allow  re-grouping  of  spectral  classes  into  composite 
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nominal  classes,  and  to  allow  manual  re-assignment  of  printer  symbols 
among  these  nominal  classes  for  more  meaningful  display  of  interpreta- 
tions. 

3.2.8  Getting  Started  with  EDITOR  and  TENEX 

In  this  section  a  brief  description  of  how  to  run  EDITOR 
under  TENEX  will  be  presented.  TENEX  is  a  time-sharing  system  running 
on  a  modified  DEC  SYSTEM-! 0  as  developed  by  Bolt,  Beranek  and  Newman 
(BBN).  [19]  When  connecting  to  a  TENEX  system,  the  user  should  estab- 
lish a  character-at-a-time,  full-duplex  connection.  Commands  to  TENEX 
may  be  abbreviated  and  if  followed  by  ESC  (ESCAPE  or  ALT  MODE  on  some 
terminals)  the  remainder  of  the  command  will  be  typed  out  followed  by 
a  prompt  for  any  parameters  required.  All  commands  should  be  terminated 
by  a  carriage  return. 

To  begin  an  EDITOR  session,  one  must  first  LOGIN.  The  se- 
quence is  LOG  <user  code>  <password>  <account  number>  followed  by  a 
carriage  return.  (The  password  will  not  be  echoed  to  preserve  secur- 
ity.) If  a  space  ic.  substituted  for  <account  number>,  that  account  num- 
ber associated  with  the  user  code  given  will  be  assumed  by  default. 

To  get  an  ERTS  image  tape  mounted,  a  request  must  be  made  to 
the  operator  who  should  be  logged  in  as  OPERATOR.  To  do  this,  the  LINK 
command  is  used.  The  LINK  command  links  two  terminals  together  so  that 
whatever  is  typed  on  either  may  be  viewed  on  both.  The  syntax  of  the 
command  is  LINK  OPERATOR.  To  communicate  with  the  operator,  a  message 
should  then  be  typed.  Each  line  of  this  message  must  begin  with  a  semi- 
colon (";")  to  indicate  to  TENEX  that  the  line  is  not  a  command.  The 
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request  should  be  to  mount  a  specific  tape,  without  write  ring  (to  in- 
sure that  the  data  on  the  tape  cannot  be  accidentally  over-written),  on 
a  9-track  or  7-track  drive.  Once  the  operator  has  mounted  the  tape,  he 
will  inform  the  user  of  the  drive.  The  tape  drives  are  MTA0:,  MTA1 : , 
etc.  The  drive  used  should  be  remembered,  since  it  will  be  needed  later 
in  using  EDITOR  to  read  the  tape.  When  discussion  with  the  operator  has 
been  completed,  the  BREAK  command  should  be  used  to  break  the  link. 

In  order  to  secure  use  of  the  tape,  two  more  TENEX  commands 
must  be  used.  The  ASSIGN  <tape  drive>  (e.g.,  ASSIGN  MTA0:)  command 
assigns  the  tape  drive  to  the  user's  job,  preventing  other  users  from 
accessing  it.  The  MOUNT  <tape  drive>  (e.g.,  MOUNT  MTA0:)  positions  the 
tape  for  reading. 

To  evoke  EDITOR,  the  user  simply  types  <RAY>  EDITOR  followed 
by  one  or  two  carriage  returns.  The  two-carriage-returns  form  provides 
a  shorter  introductory  herdld  and  is  recommended  for  regular  users.  The 
user  is  then  at  the  top  command  level  of  EDITOR  and  may  enter  commands 
as  detailed  in  the  following  section.  A  list  of  legitimate  commands  is 
available  by  entering  question  mark  ("?")  followed  by  a  carriage  return 
(CR).  At  the  top  level,  any  command  may  be  abbreviated—as  short  as 
required  to  distinguish  it  from  all  other  commands—and  followed  by  a 
carriage  return  or  ESC.  If  followed  by  ESC,  the  remainder  of  the  abbre- 
viated command  will  be  completed  by  EDITOR.  In  certain  of  the  analysis 
modules  within  EDITOR,  a  somewhat  different  command  structure  is  used  in 
which  the  user  is  prompted  for  short  inputs  by  a  series  of  questions. 
To  return  to  the  top  level  of  EDITOR  from  any  lower-level  module,  an 
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exclamation  point  ("!"  and  CR)  is  used  instead  of  a  command.  (One  is 
at  the  top  level  of  EDITOR  when  one  is  prompted  by  a  single  exclama- 
tion point. ) 

Most  of  the  EDITOR  modules  require  TENEX  disk  files  for  input 
and/or  output  of  data.  The  simplest  form  of  a  file  designator  is  FILE- 
NAME.EXTENSION  where  FILENAME  and  EXTENSION  are  identifiers  and  the 
period  is  part  of  the  file  designator.  Note  that  the  period  (".")  is 
always  required  between  FILENAME  and  EXTENSION. 

In  entering  commands  through  a  terminal,  errors  will  occa- 
sionally be  made.  To  delete  a  single  character,  CONTROL-A  is  used  and 
to  delete  an  entire  line,  CONTROL-X  is  used.  Similar  conventions  are 
also  available  at  the  level  of  TENEX  commands.  For  example,  TENEX  pro- 
vides CONTROL-C  for  terminating  any  program  and  returning  to  TENEX 
EXEC.  (In  some  cases,  two  CONTROL-C  characters  must  be  entered.)  It 
should  be  noted  that,  wherever  one  might  be  in  EDITOR,  the  CONTROL-C 
completely  exits  EDITOR.  However,  this  is  the  only  way  to  terminate  a 
long  program,  such  as  a  large  cluster  or  classify,  which  was  somehow 
initiated  erroneously. 

Finally,  when  the  EDITOR  session  is  complete,  EDITOR  is  exited 
using  the  QUIT  EDITOR  command.  If  a  tape  was  used,  the  user  should 
relinquish  control  of  the  tape  drive.  First  the  DEASSIGN  <tape  drive> 
(e.g.,  DEASSIGN  MTA0:)  command  is  used  to  release  exclusive  control  of 
the  tape  drive.  Next  the  UNLOAD  <tape  drive>  (e.g.,  UNLOAD  MTA0:) 
command  is  used  to  rewind  the  tape  and  position  it  for  dismounting,  hence 
signaling  the  operator  that  he  may  remove  the  tape  and  put  it  away. 
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Thus,  it  is  not  necessary  to  LINK  to  the  operator  again  at  the  end  of  an 
EDITOR  session. 

Before  logging  out,  the  user  should  check  to  see  if  there  are 
any  files  on  TENEX  disk  that  will  be  no  longer  needed.  A  listing  of  all 
files  in  the  user's  directory  may  be  obtained  via  the  TENEX  command 
DIRECTORY.  It  is  prudent  to  delete  files  since  all  disk  space  used  costs 
money,  and  since  once  a  directory  is  filled  to  its  allocation  (the 
DSKSTAT  command  will  indicate  the  space  remaining)  no  new  files  may  be 
created. 

To  delete  files,  the  command  DELETE  <file  name>  is  used. 
Where  it  is  desirable  to  delete  all  files  of  a  common  FILENAME  (or  prefix 
within  the  FILENAME. EXTENSION  convention),  this  may  be  done  conveniently 
by  using  the  DELETE  FILENAME.*  form  of  the  DELETE  command.  This  option 
for  TENEX  file  deletion  should  be  noted  by  EDITOR  users  since,  if  care 
is  given  to  naming  files  within  EDITOR  sessions,  file  directory  mainte- 
nance can  be  greatly  simplified.  In  any  case,  if  obsolete  TENEX  files 
are   being  deleted  in  order  to  regain  disk  space  needed  for  immediate 
use,  the  TENEX  EXPUNGE  command  should  be  issued  by  the  user  following 
all  DELETES. 

When  the  TENEX  session  is  complete,  the  LOGOUT  command  is  used 
to  exit  the  system  before  closing  the  connection.  At  this  point,  the 
TENEX  system  will  report  to  the  user  the  total  amount  of  processing  (CPU) 
time  and  terminal-connect  time  used  during  the  complete  session. 
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3.2.9  Summary  of  EDITOR  Commands 

In  this  section,  a  brief  description  of  the  various  EDITOR 
commands  and  the  modules  they  evoke  will  be  given.  The  use  of  most  of 
these  commands  is  illustrated  in  the  example  EDITOR  sessions  given  in 
the  Appendix.  Several  commands  have  already  been  described  partially 
in  Section  3.2.7  above.  The  commands  are  discussed  here  in  alphabetical 
order  to  assist  subsequent  referencing. 

3.2.9.1  CLASSIFY 

The  CLASSIFY  command  is  used  to  perform  statistical  classifica- 
tion of  an  image  data  file  using  a  TENEX  routine.  The  number  of  channels 
will  be  requested.  At  present,  only  4-  or  8-channel  data  may  be  handled 
so  "4"  or  "8"  should  be  entered.  The  user  will  be  prompted  for  other 
required  inputs.  For  large-scale  classifications,  the  ILLIAC  IV  may  be 
used  via  the  EXECUTE  ANALYSIS  ON  ILLIAC  IV  command  described  below. 

3.2.9.2  CLIP  A  WINDOW 

The  CLIP  A  WINDOW  command  is  used  to  create  subwindows  of 
windows  in  an  image  file  already  on  TENEX  disk.  A  question  mark  ("?") 
upon  entering  CLIP  lists  the  commands.  CLIP  may  be  used  for  editing  4- 
or  8-channel  raw  or  categorized  data  files.  The  subwindows  to  be  created 
may  be  specified  in  terms  of  new  lines  and  columns  (contained  within  the 
input  window)  and/or  additional  image  sampling.  Commands  within  the 
CLIP  A  WINDOW  module  allow  selection  of  the  image  window  to  be  edited, 
and  the  new  coordinates  and  sampling  increments  to  be  used  in  creating 
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the  new  subwindows.  After  all  window  editing  parameters  have  been 
specified,  a  command  is  available  to  write  the  output  file  and  return 
to  the  top  level  of  EDITOR. 

3.2.9.3  CLUSTER 

The  CLUSTER  command  is  used  to  perform  cluster  analysis  on 
image  data  using  TENEX  processing.  The  number  of  channels  will  be 
requested.  At  present,  only  4-  or  8-channel  data  may  be  handled  so 
"4"  or  "8"  must  be  entered.  The  user  will  be  prompted  for  other  required 
inputs.  For  large  cluster  analyses,  the  ILLIAC  IV  may  be  used  via  the 
EXECUTE  ANALYSIS  ON  ILLIAC  IV  command  (again,  described  below). 

3.2.9.4  ENTER  STATISTICS  FILE  EDITOR 

The  ENTER  STATISTICS  FILE  EDITOR  command  enters  the  statistics 
editor  module.  Subcommands  are  of  the  form  one  letter  followed  by  a 
carriage  return.  A  list  of  commands  may  be  displayed  by  typing  question 
mark  ("?")  followed  by  carriage  return.  Within  this  module,  the  summary 
statistics  for  sets  of  spectral  categories  generated  by  cluster  analyses 
may  be  displayed.  Also,  the  statistics  of  individual  categories  of 
various  files  may  be  grouped  in  new  combinations  to  form  new  statistics 
files.  In  addition,  the  statistics  of  several  categories  may  be  pooled 
together  to  create  a  composite  spectral  "signature."  Finally,  where  a 
classification  is  to  be  done  using  the  ILLIAC  IV,  the  variance-covariance 
matrices  of  the  appropriate  statistics  file  may  be  inverted. 
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3.2.9.5  EXECUTE  ANALYSIS  ON  ILLIAC  IV 

The  EXECUTE  ANALYSIS  ON  ILLIAC  IV  command  enables  automatic 
submission  and  retrieval  of  ILLIAC  IV  batch  jobs  from  a  remote  TENEX 
site.  The  ARPA  Network  is  used  to  submit  the  jobs  and  retrieve  the 
results  from  the  ILLIAC  IV  at  NASA-Ames.  To  use  the  ILLIAC  IV,  the  user 
must  have  a  valid  user  code  and  password  at  I4-TENEX  (the  TENEX  front- 
end  subsystem  of  the  ILLIAC  IV  complex).  EDITOR  will  ask  for  these  be- 
fore allowing  any  commands  to  be  entered.  At  present,  both  cluster 
analysis  and  statistical  classification  procedures  may  be  performed  on 
the  ILLIAC  IV  remotely  through  EDITOR. 

3.2.9.6  IDENTIFY  A  WINDOW 

The  IDENTIFY  A  WINDOW  command  displays  information  extracted 
from  the  header  of  an  image  file.  This  information  includes  the  image 
file  type,  the  line  and  column  coordinates  of  all  windows  in  the  file, 
the  line  and  column  sampling  increments  used  in  creating  all  windows  in 
the  file,  the  number  of  channels  of  the  image  data,  the  cluster  analysis 
parameters  (DELTA  or  separability  threshold  and  the  NUMBER  OF  CLUSTERS), 
and  certain  descriptive  information  supplied  by  NASA  for  the  ERTS  image 
from  which  the  data  have  been  extracted. 

3.2.9.7  INDICATE  A  PRIORI  PROBABILITIES 

The  INDICATE  A  PRIORI  PROBABILITIES  command  allows  creation  of 
a  file  of  a  priori  probabilities  of  various  categories  for  classification. 
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3.2.9.8  MODIFY  WINDOW  HEADER 

The  MODIFY  WINDOW  HEADER  command  allows  insertion  of  neces- 
sary parameters  (NUMBER  OF  CLUSTERS  and  DELTA  or  separability  threshold) 
to  be  used  for  a  particular  cluster  analysis.  Also  there  the  alpha- 
numeric information  recorded  in  the  header  of  an  image  file  may  be  modi- 
fied. Once  all  desired  changes  are  made,  the  CLOSE  FILE  command  should 
be  issued  to  return  to  the  top  level  of  EDITOR. 

3.2.9.9  PRINT  A  WINDOW 

The  PRINT  A  WINDOW  command  allows  display  of  image  windows. 
The  display  may  be  directed  either  to  the  user's  terminal  or  to  a  back-up 
disk  file  for  later  printing  using  a  conventional  line  printer.  Raw 
data  windows  are  displayed  with  an  overprinting  routine  allowing  up  to 
eleven  (11)  levels  of  gray-scale.  A  histogram  of  the  distribution  of 
spectral  intensities  for  any  spectral  channel  may  also  be  displayed.  For 
categorized  windows,  the  character  a-ssigned  to  each  class  may  be  selected 
manually  as  any  printing  character  to  achieve  more  meaningful  displays. 
This  feature  facilitates  the  re-grouping  of  spectral  classes  by  nominal 
categories  for  tabulation  and  display  purposes. 

NOTE:  When  PRINT  A  WINDOW  output  is  directed  to  a  terminal, 
the  line  width  used  for  displays  of  windows  is  the  line  width  assigned 
to  that  terminal  by  TENEX.  The  default  value  is  72  characters,  but  this 
may  be  changed  by  the  TENEX  WIDTH  <number>  command,  where  <number>  is 
the  width  desired.  Since  WIDTH  is  a  TENEX  command  it  must  be  performed 
before  entering  EDITOR. 
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3.2.9.10  QUIT  EDITOR 

The  QUIT  EDITOR  command  provides  for  a  proper  exit  from  the 
top  EDITOR  command  level  back  to  TENEX  EXEC.  The  QUIT  command  should 
always  be  used  whenever  possible  (instead  of  CONTROL-C)  to  exit  EDITOR. 

3.2.9.11  RETRIEVE  WINDOW  FROM  TAPE 

The  RETRIEVE  WINDOW  FROM  TAPE  command  is  used  to  copy  windows 
from  a  specific  tape  (that  has  already  been  mounted)  into  a  disk  file. 
Where  desired,  line  and  column  sampling  increments  should  be  supplied 
via  the  SAMPLE  command.  Then  the  coordinates  of  the  windows  to  be 
retrieved  are  entered  using  the  INSERT  COORDINATES  command.  Additional 
facilities  are  available  to  allow  the  user  to  delete  window  coordinates 
entered  in  error. 

Prior  to  any  EDITOR  analysis  of  a  specific  ERTS  image,  specific 
geographic  areas  of  interest  are  known  to  the  user,  typically,  only  in 
terms  of  latitude  and  longitude  boundaries.  On  the  other  hand,  all  image 
windows  to  be  retrieved  from  tape  by  EDITOR  must  be  specified  in  terms 
of  image  line  and  column  coordinates.  Therefore  as  a  convenience  to  the 
user,  there  is  included  within  the  RETRIEVE  module  a  subprocedure  for 
transformation  of  window  corners  specified  in  terms  of  latitude  and 
longitude  coordinates  into  "nearest-approximation"  image  pixels  identi- 
fied by  image  line  and  column  numbers.   (The  inverse  transformation  is 
also  possible;  that  is,  for  image  pixels  identified  by  line  and  column 
coordinates,  approximate  latitude  and  longitude  coordinates  may  be  com- 
puted. ) 
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4.   PATTERN  RECOGNITION:  A  BASIS  FOR 
REMOTE  SENSING  DATA  ANALYSIS 


4.1  Classification  Procedure 


Pattern  recognition  plays  a  central  role  in  numerically 
oriented  remote  sensory  systems.   It  provides  an  automatic  procedure 
for  deciding  to  which  class  any  given  ground  resolution  elements 
should  be  assigned. 

Figure  6  shows  a  model  of  a  general  pattern  recognition 
system.   In  the  LARS  [7]  context  the  receptor  or  sensor  is  usually  a 
multispectral  scanner.  For  each  ground  resolution  element  the  sensor 
produces  n_  numbers  or  measurements  corresponding  to  the  n   channels  of 
the  scanner.   It  is  convenient  to  think  of  the  n   measurements  as  defin- 
ing a  point  in  n-dimensional  Eucl idean  Space  which  is  referred  to  as 
the  Measurement  Space.  Any  particular  measurement  can  be  represented  by 
the  vector: 

~xl 


X  = 


The  feature  extractor  transforms  the  n-dimensional  measurement  vector 

i 
into  an  n  -dimensional  feature  vector.  The  decision  maker  in  Figure  7 

performs  calculations  on  the  feature  vectors  presented  to  it  and,  based 
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Figure  6  A  Simplified  Model  of  a  Pattern  Recognition  Syst 
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Figure  7  Decision  Regions  and  Surfaces 
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Figure  8  A  Pattern  Classifier  Defined  in  Terms  of  Discriminant  Functi 
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upon  a  decision  rule,  assigns  the  "unknown"  data  point  to  a  particular 
class. 

4.1.1  Discriminant  Functions:  Quantifying  the  Decision  Procedure 

Patterns  arising  in  remote  sending  problems  exhibit  some  ran- 
domness due  to  the  randomness  of  nature.  Vectors  from  the  same  class 
tend  to  form  a  "cloud  of  points"  as  shown  in  Figure  7.  The  job  of  the 
pattern  classifier  is  to  divide  the  feature  space  into  decision  regions, 
each  region  corresponding  to  a  specific  class.  Any  data  point  falling 
in  a  particular  region  is  assigned  to  the  class  associated  with  that 
region.  The  surface  separating  the  decision  regions  are  known  as  deci- 
sion surface.  Designing  a  pattern  recognizer  really  boils  down  to  form- 
ing a  procedure  for  determining  the  decision  surface  so  as  to  optimize 
some  performance  criterion,  such  as  maximizing  the  frequency  of  correct 
classification. 

These  concepts  can  be  put  on  a  quantitative  basis  by  intro- 
ducing Discriminant  Functions.  Assume  there  are  m  pattern  classes.  Let: 

g1  (x),  g2  (x),...,  gm  (x) 

be  scaler  single  valued  functions  of  x  such  that 

gi  (x)  >  g,  (x)  for  all  x 

in  the  corresponding  region  corresponding  to  the  i   class  (J  f   i).  If 
the  discriminant  functions  are  continuous  across  the  decision  boundaries, 
the  decision  surface  are  given  by  equation  of  the  form. 


9^   (x)  -  g.   (x)  =  0 
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A  pattern  classifier  can  then  be  represented  by  the  block  diagram  of 
Figure  8.    By  taking  this  approach  the  pattern  classifier  design 
problem  is  reduced  to  the  problem  of  how  to  reflect  the  discriminant 
functions  in  an  optimal  fashion. 

4.1.2  "Training"  the  Classifier 

In  some  cases  it  is  possible  to  select  discriminant  functions 
on  the  basis  of  theoretical  considerations  or  experience.  More  com- 
monly the  discriminant  functions  are  based  upon  a  set  of  training  pat- 
terns. Training  patterns  which  are  typical  of  those  to  be  classified 
are  "shown  to  the  classifier  together  with  the  identity  of  each  pattern 
and  based  on  this  information  the  classifier  establishes  its  discriminant 
functions. 

gi  (x),  i  =  1,2,  ....  n. 

Example:   Consider  a  two  dimensional,  two  class  problem  in  which  the 
discriminant  functions  are  assumed  to  have  the  form 

g1  (x)  =  an  x]  +  9]?   x2  +  b] 

g2  (x)  =  a21  x1  +  g22  x2  +  b2 

Then  g,  (x)  -  g2  (x)  =  0  is  the  equation  of  a  straight  line  dividing  x, , 
x«  plane.  Given  a  set  of  training  patterns,  how  should  the  constants 
a,-.,  a,2»  b-i  >  etc.  be  chosen?  It  can  be  proven  that  if  the  training  pat- 
tern are  indeed  separable  by  a  straight  ine,  then  the  following  procedure 
will  converge.  [8] 
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Initially  select  a's  and  b's  arbitrarily.  For  example,  let: 
all  =  a12  =  bl  =  ] 

Oiry-t        ~       ®00       ~  9   _   ~ 

Then  take  the  first  training  pattern  (say  it  is  from  w, ,  i.e.,  from 
class  1)  and  calculate  g,  (x)  and  g2  (x).  If  g,  (x)  >  g2  (x)  the  deci- 
sion is  correct;  go  on  to  the  next  training  sample.  If  g-,  (x)  <  g2  (x) 
a  wrong  decision  would  be  made.  In  this  case,  alter  the  coefficients  so 
as  to  increase  the  discriminant  function  associated  with  the  correct 
class  and  decrease  the  discriminant  function  associated  with  the  incor- 
rect class.  If  x  is  from  to-,  but  co2  was  decided,  let 

i  i 

a -in  -  aii   cxXi       ^9i  —  ^9i  ™  o'X  i 

i  i 

3i2  =  3-j  n  ■  CXX  o  ^99  =  ^99  ~  CXX^ 

i  i 

b-,  =  b,  +  a         bp  =  b0  -  a 

where  a  is  a  convenient  positive  constant.  If  x  is  from  w^  but  w-,  was 
decided,  change  the  signs  in  the  above  equation  so  as  to  increase  g2  and 
decrease  g-,.  Then  go  to  the  next  training  pattern.  Cycle  through  the 
training  patterns  until  all  are  correctly  classified. 

4.1.3   The  Statistical  Approach 

Remote  sensing  is  typical  of  many  practical  applications  of 
pattern  recognition  for  which  statistical  methods  are  appropriate  in  the 
following  respects: 

1.  The  data  exhibit  many  incidental  variations  (noise)  which 
tend  to  obscure  differences  between  the  pattern  classes. 
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2.  The  pattern  classes  of  interest  may  actually  overlap  in 
the  measurement  space  (may  not  always  be  discriminable),  suggesting  the 
use  of  an  approach  which  leads  to  decisions  which  are  "most  likely" 
correct. 

Statistical  pattern  recognition  techniques  often  make  use  of 
the  probability  density  functions  associated  with  the  pattern  classes 
(including  the  approach  to  be  described  here).  However,  the  density 
functions  are  usually  unknown  and  must  be  estimated  from  a  set  of  train- 
ing patterns.  In  some  cases,  the  form  of  the  density  functions  is 
assumed  and  only  certain  parameters  associated  with  the  functions  are 
estimated.  Such  methods  are  called  "parametric. "   Methods  for  which 
not  even  the  form  of  the  density  functions  is  assumed  are  called  "non- 
parametric.  "   The  parametric  case  requires  more  a  priori  knowledge  or 
some  basic  assumptions  regarding  the  nature  of  the  patterns.  The  non- 
parametric  case  requires  less  initial  knowledge  and  fewer  assumptions 
but  is  generally  more  difficult  to  implement.  Let  there  be  m  classes 
characterized  by  the  conditional  probability  density  functions. 

p  (X^.)     i  =  1,2,.  .  .,  m  (1) 

The  function  p  (XJoo.)  gives  the  probability  of  occurrence  of  pattern  X, 
given  that  X  is  in  fact  from  class  i. 

An  important  assumption  in  the  LARSYS  (Purdue  University) 
algorithms  is  that  the  p  (X|gj.  )  are  each  multivanate  Gaussian  (or  normal) 
distributions.  This  is  a  parametric  assumption  which  leads  to  a  form  of 
classifier  which  is  relatively  simple  to  implement.  Under  this  assumption, 
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a  mean  vector  and  covariance  matrix  are  sufficient  to  characterize  the 
probability  distribution  of  any  pattern  class. 

Returning  to  the  problem  of  how  to  specify  the  discriminant 
functions,  an  approach  based  on  statistical  decision  theory  is  taken. 
A  set  of  less  functions  is  defined: 

X   (i|j)  i,j  =  1,2,.  .  .,  m  (2) 

where  A  (i|j)  is  the  loss  (or  cost)  incurred  when  a  pattern  is  classified 
into  class  j_  when  it  is  actually  from  class  j_.  If  the  pattern  classifier 
is  designed  so  as  to  minimize  the  average  (expected)  loss,  then  the 
classifier  is  said  to  be  Bayes  optimal .  For  a  given  pattern  X,  the 
expected  loss  resulting  from  the  decision  Xeo).  is  given  by 

m 

Lx  (D  =  I   A  (i|j)  pU  |x)  (3) 

j  =  l  J 

where  p(w.|x)  is  the  probability  that  a  pattern  X  is  from  class  j. 
Applying  Bayes'  rule,  i.e., 

p  (x,Wj)  =  p  (x|w.)  P  (wj)  =  p(o).jx)  p  (x)  (4) 

the  expected  loss  can  be: 

Lx  (D  =  I     A  Ci  id)  P(x|a).)  p(u).)/p(x)  (5) 

j  =  l  J     J 

where  p  (go.)  is  the  priori  probability  of  to., 
j  j 

Note  the  minimizing  Lx  (i)  with  respect  to  i_  is  the  same  as 
maximizing  -  Lx  (i).  Thus  a  suitable  set  of  discriminant  functions  is 

9,-  (x)  =  -Lx  (1)   i  =  1,2,.  .  .,m  (6) 
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A  simple  (a  reasonable)  loss  function  is 
X  (i|j)  =  0         i  =  j 


A  (i|j)  =  1         i  t   J 

m 
Then     g.(x)  =  -  \     p(x|u>.)  p(ui.)/p(x)  (7) 

J7i 
Note  that  from  any  set  of  discriminant  functions,  another  set  of  dis- 
criminant functions  can  be  formed  by  taking  the  same  monotonic  functions 
of  each  of  the  original  discriminant  functions.  For  example,  if: 

gi  (x),  1  =  1,2,  •  .  -,m 

is  a  set  of  discriminant  functions,  then  so  are  the  sets 

9,-  U)  =  9,-  (x)  +  constant  i  =  1,2,.  .  .  ,m 

and 

9-|  (x)  =  log  [g^x)]      i  =  1,2,.  .  .  ,m 

Examing  (7)  note  that  p  (x)  is  not  a  function  of  i  so  it  is  just  as  well 

to  maximize. 

m 

9,-  (x)  =  -  I         P  (xLo)  p  (o>.)  =  -  [p(x)  -  p(x  to. )  p  (u>.)] 
i        j=1       I  J      J  «  l      i 

But  this  is  maximum  if 

q        (x)  =  p  (x|ooi )  p  (w.) 

is  maximum.  Thus,  the  decision  rule  is: 
decide   x  ecu.  if  and  only  if 


p  (x|uji )  p  (x^  >  p  (x|oa.)  p  (go.)  for  all  j 
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that  is  deciding  x  eu.  if  g.  (x)  =  g.  (x)  and  i  >  j.  This  is  commonly 
referred  to  as  the  Maximum  Likelihood  Decision  Rule. 

Consider  the  maximum  likelihood  discriminant  function  as  it 
applies  to  remote  sensing.  The  p  (w. )  represent  the  priori  probability 

t  h 

of  the  i   class.  This  can  often  be  estimated.  Taking  agricultural 
crop  types  as  an  example,  the  p  (to. )  may  be  estimated  from  previous  year 
yields,  seed  sales  records  or  statistical  reporting  service  information. 
The  densities  p  (xju).),  on  the  other  hand,  generally  have  to  be  esti- 
mated from  training  samples. 

The  assumptions  upon  which  the  classification  algorithms  are 
based  is  that  p  (x|w. )  is  a  multivariate  Gaussian  probability  density 
function.  Examining  the  maximum  likelihood  decision  criterion  in 
one-dimensional  Gaussian  case  will  serve  both  as  a  review  of  Gaussian 
density  functions  and  as  a  means  of  illustrating  the  principles  of  pat- 


tern classification.  In  this  case  (e.g.,  one  spectral  channel) 
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when  u.  =  E  [x]  and  a.  =  E  [(x-u^)  ]  are  the  mean  and  variance  for 

2 
class  i.  In  practice  u.  and  a.  are   unknown  and  must  be  estimated  from 

training  samples.  From  statistical  theory,  we  have 

n, 
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I       (x.-m.)' 
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(n.  =  number  of  training  patterns  in  class  1)  are  unbiased  estimators  of 
the  mean  and  variance.  Thus,  the  estimated  density  function  is: 

2 


P  (xjoj^  =  1 


1/2 
(2tt)   Si 


exp 


-1/2  <X-mi>' 


1      _J 

Following  the  decision  theory  approach  the  discriminant  function  is 


9i  (x) 


p(u).j) 


exp 


(2ir) 


1/2 


Si 


-  1/2  (x-mi}< 


S2 
bi 


and  since  a  monotonic  function  of  a  discriminant  function  may  also  be 
used  as  a  discriminant  function,  we  shall  take  the  logarithm  of  the 
previous  function  to  obtain: 

(  \2 

g.  (x)  =  log  p  (u>.)  -  1/2  log  2  tt  -  logs.-  1/2  u  "  V 


g.'    (x)  =  log  p  (a).)  -  log  Si  -  1/2  ^"m1)' 


Thus  the  decision  rule  becomes: 


Decide  Xcw.  if  and  only  if 


2  2 

(x-m.)  (x-m.) 

log  p(o)i)  -log  S-i/2  — ^ —  >  log  p(u.)  -  log  S-.   1/2  —^ — 

si  S. 


In  the  two  dimensional  case 
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and 
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where 


yil  =  E  l-xl  |0Ji-''  yi2  =  E  ^x2  |a)i-' 
Qijk  =  E[(xi  "  »*1J>  (xk  -  yik>h] 


by  defining  a  mean  vector  and  covariance  matrix 


j,k  =  1,2 

i  =  1,2,  .  .  .R 


Ui  = 


II- 


Mil 
yi2 


°ill  ail2 


ai21  ai22 


The  density  can  be  rewritten  in  the  simple  form: 


P  (xjo)^  =  1_ 


2ir|I.| 


1/2 


exp 


1/2  (x 


U,)  T  I-1  (x  -Uj) 
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Wh 


ere  |J.|  is  the  determinant  of  J.  and  (x-U.)  is  the  transpose  of 


(x  -  U.).  Thus  this  matrix  formulation  holds  for  n_  dimensions  as  well 
as  for  two  dimensions. 

4.2  Clustering  Procedure 

Clustering  is  a  data  analysis  technique  by  which  one  attempts 
to  determine  the  "natural"  relationships  in  a  set  of  observations  or 
data  points.  It  is  sometimes  referred  to  as  unsupervised  classifica- 
tion because  the  end  product  is  generally  a  classification  of  each 
observation  into  a  "class"  which  has  been  established  by  the  analysis 
procedure,  based  on  the  data,  rather  than  by  the  person  interested  in 
the  analysis.  Clustering  has  been  applied  as  a  means  of  data  compres- 
sion and  for  the  purpose  of  determining  differentiating  characteristics 
in  complex  data  sets.  An  increasingly  important  application  is  unsuper- 
vised classification,  in  which  the  clustering  algorithm  determines  the 
classes  based  on  the  clustering  tendencies  in  the  data.  Essentially,  the 
definition  of  a  clustering  algorithm  depends  on  the  specification  of  two 
distance  measures:  a  measure  of  distance  between  data  points  or  indi- 
vidual observations,  and  a  measure  of  a  distance  between  groups  of  obser- 
vations. Figure  9  is  a  block  diagram  for  a  typical  clustering  algorithm. 
The  point-to-point  distance  measure  is  used  in  the  step  labeled  "assign 
each  vector  to  nearest  cluster  center."  The  distance  between  groups  of 
points  (clusters,  in  this  case)  is  calculated  in  the  step  "compute 
separability  information."  Euclidean  distance,  the  most  familiar  point- 
to-point  distance  measure,  is  defined  for  two-n-dimensional  points  or 


61 


Initialize 

cluster 

centers 


Assign  each  vector 
to  nearest 
cluster  center 


Calculate  means  of 

new  "clusters" 

(new  cluster  centers) 


Yes 
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Figure  9  Clustering  Algorithm 


M  ? 

Euclidean  distance:  D  =  [£  (x.-y.)  ] 

i=l   1  n 


62 
1/2 


Several  alternatives  are  available  as  candidate  measures  of  distance 
between  clusters.  One  possibility  is  the  divergence  or  transformed 
divergence  used  for  feature  selection.  In  LARSYS,  a  measure  called 
"Swain-Fu  distance"  has  been  implemented,  which  compares  the  separation 
of  cluster  centers  to  the  dispersion  of  the  data  in  the  clusters.  The 
dispersion  of  the  data  in  a  cluster  is  measured  in  terms  of  the 
"ellipsoid  of  concentration"  associated  with  the  cluster. 

Ellipsoid  of  concentration:  Let  the  random  vector  X  have  a 
distribution  with  mean  vector  U  and  covariance  matrix  I  =   [<?■,•■;]•  If 
Z  is  another  random  vector  uniformly  distributed  over  the  volume  of  the 
ellipsoid  given  by 

Q  (Z)  =  I       I       \hA     (Z.-U.)  (Z.-U.)  =  n  +  2 
1-1  J-l   m  J  3 

where  n  is  the  number  of  components  in  1,   and  |T.  .  |  is  the  cofactor  of 
ff. .,  then  Z  also  has  zero  mean  and  covariance  matrix  Z.  The  ellipsoid 
Q   is  called  the  ellipsoid  of  concentration  of  the  distribution  of  X. 

Consider  two  clusters  and  their  respective  ellipsoids  of  con- 
centration as  shown  in  Figure  10.  D-.^  is  the  distance  between  the 
cluster  centers.  D-.  is  the  distance  from  the  center  of  cluster  1  to  the 
surface  of  its  ellipsoid  of  concentration  along  the  line  connecting  the 
cluster  centers  similarly  to  D«.  In  terms  of  these  distances,  D, ,  Dp, 
D-|2»  the  Swain  in  distance  is  proven  by 

A=   D12 


Dl  +D2 
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In  terms  of  the  cluster  centered  (cluster  means)  and  the  covariance 
matrices  associated  with  the  clusters,  the  Swain-Fu  distance  can  be 
expressed  as 


A  =  TzTcZ 


/c7  +  /c\. 


where    Ck  =  tr  {£~]  (Uj-u2)  (uru2)T} 

tr  {A}  =  trace  of  matrix  A 

J,  =  covariance  matrix  for  cluster  K 

U.  =  mean  vector  for  cluster  K. 

An  illustration  of  how  the  clustering  algorithm  works  refers 
to  Figure  10.  The  first  step  is  to  select  initial  cluster  centers.  The 
analyst  must  specify  how  many  clusters  are  to  be  isolated;  the  algorithm 
determines  (arbitrarily)  where  the  initial  centers  are  to  be  located  (the 
final  results  are  relatively  insensitive  to  the  initial  selection.  Each 
data  point  is  then  labeled  as  "belonging"  to  the  nearest  cluster  center 
(using  Euclidean  distance),  effectively  creating  a  cluster  of  data  points 
associated  with  each  center.  The  boundaries  between  clusters  are  formed 
by  the  lines  (planes  in  n-dimensional  space)  which  are  the  perpendicular 
bisectors  of  the  lines  connecting  the  centers.  Next,  new  cluster  cen- 
ters are  calculated.  The  new  center  for  each  cluster  is  the  mean  vector 
of  all  points  just  assigned  to  that  cluster.  A  check  to  see  whether  the 
algorithm  has  achieved  the  final  results,  which  is  the  care  when  the  new 
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Figure  10  Separability  of  Clusters 
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cluster  centers  are  identical  with  the  previous  centers.  When  no 
further  change  is  detected,  the  pair  wise  distances  (Swain-Fu  distance) 
between  the  resulting  clusters  are  computed  and  all  results  are  printed 
for  evaluation  by  the  analyst. 
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